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ABSTRACT 

Pseudo-transient  continuation  (^tc)  is  a  well-known  and  physically  motivated  tech¬ 
nique  for  computation  of  steady-state  solutions  of  time-dependent  partial  differential 
equations.  Standard  globalization  strategies  such  as  line  search  or  trust  region  meth¬ 
ods  often  stagnate  at  local  minima,  ^tc  succeeds  in  many  of  these  cases  by  taking 
advantage  of  the  underlying  PDE  structure  of  the  problem.  Though  widely  employed, 
the  convergence  of  ^tc  is  rarely  discussed.  In  this  paper  we  prove  convergence  for  a 
generic  form  of  ^tc  and  illustrate  it  with  two  practical  strategies. 


^  This  author’s  work  was  supported  by  NSF  grant  #DMS-9321938. 

^  This  author’s  work  was  supported  in  part  by  NSF  grant  ECS-9527169,  and  by  NASA  contracts 
NAGI-1692  and  NASl-19480,  the  latter  while  the  author  was  in  residence  at  the  Institute  for  Com¬ 
puter  Applications  in  Science  and  Engineering  (ICASE),  NASA  Langley  Research  Center,  Hampton, 
VA  2.3681-0001. 


i 


1.  Introduction.  Pseudo-transient  continuation  (^tc)  is  a  method  for  compu¬ 
tation  of  steady-state  solutions  of  partial  differential  equations  framed  in  a  time- 
dependent  setting,  x'  =  F{x).  Here 


and  spatial  variables  and  derivatives  are  contained  in  the  nonlinear  term  F{x).  March¬ 
ing  out  the  steady  state  by  following  the  physical  transient  may  be  unnecessarily  time 
consuming  if  the  intermediate  states  are  not  of  interest.  On  the  other  hand,  Newton’s 
method  for  F(a:)  =  0  alone  will  usually  not  suffice,  as  initial  iterates  sufficiently  near 
the  root  are  usually  not  available.  Standard  globalization  strategies  [12,  17,  24],  such 
as  line  search  or  trust  region  methods  often  stagnate  at  local  minima  of  ||F||  [19]. 
This  is  particularly  the  case  when  the  solution  has  complex  features  such  as  shocks 
that  are  not  present  in  the  initial  iterate  (see  [23],  for  example),  ^tc  succeeds  in 
many  of  these  cases  by  taking  advantage  of  the  PDE  structure  of  the  problem. 

1.1.  The  Basic  Algorithm.  In  the  simple  form  considered  in  this  paper,  'I'tc 
numerically  integrates  the  initial  value  problem 

(1.1)  x'  =  —V~^F{x),  a:(0)  =  xq 

to  steady  state  using  a  variable  timestep  scheme  that  attempts  to  increase  the  timestep 
as  F{x)  approaches  0.  V  is  a  nonsingular  matrix  used  to  improve  the  scaling  of  the 
problem. 

We  can  define  the  '$'tc  sequence  by 

(1.2)  =Xn-  {S~^V  -t-  F’{Xn))~'^F{Xn) 

where  F'  is  the  Jacobian  (or  the  Frechet  derivative  in  the  infinite- dimensional  situa¬ 
tion).  Algorithmically, 

Algorithm  1.1. 

1.  Set  X  =  xq  and  6  =  So-  Evaluate  F{x). 

2.  While  ||F(a:)||  is  too  large. 

(a)  Solve  (<5-ip  +  F'{x))s  =  -F{x). 

(b)  Set  X  =  X  +  s. 

(c)  Evaluate  F{x). 

(d)  Update  6. 

The  linear  equation 

(1..3)  (r'y  +  F'(:c))s  = -F’(x) 

for  the  Newton  step  that  is  solved  in  step  2a  is  often  the  discretization  of  an  elliptic 
PDE.  As  such,  it  is  typically  solved  by  an  iterative  method  which  terminates  on  small 
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linear  residuals.  This  results  in  an  inexact  method  [11],  and  a  convergence  theory  for 
’^tc  must  account  for  this.  We  have  not  yet  explained  how  6  is  updated,  nor  explored 
the  reuse  of  .Jacobian  information  from  previous  iterations.  These  options  must  be 
considered  to  explain  the  performance  of  ’I'tc  as  practiced,  and  we  take  them  up  in 
§  .3.  In  order  to  most  clearly  explain  the  ideas,  however,  we  begin  our  analysis  in  §  2 
with  the  most  simple  variant  of  ^tc  possible  and  then  extend  those  results  to  the 
more  general  situation. 

As  a  method  for  integration  in  time,  ^tc  is  a  Rosenbrock  method  ([14],  p.  223) 
if  8  is  fixed.  One  may  also  think  of  this  as  a  predictor-corrector  method,  where  the 
simple  predictor  (result  from  the  previous  timestep)  and  a  Newton  corrector  are  used. 
To  see  this  consider  the  implicit  Euler  step  from  x„  with  timestep  6„, 

(1.4)  -^(^n+l)* 

In  this  formulation  2,^+1  would  be  the  root  of 

G{0==^  +  ^nV-^F{0-Xn. 

Finding  a  root  of  G  with  Newton’s  method  would  lead  to  the  iteration 

6+1  =  6  -  (/  +  ^nV-^F(6))-^(6  +  ^ni^-^F’(6)  -  Xn) 

If  we  take  ^0  =  Xn  the  first  Newton  iterate  is 

+  6nV-'^F'{Xr.))-Hr,V-^F{Xn)  =  Xr,  -  {S^^V  +  F'{Xn))-^  F{Xr.). 

This  leaves  open  the  possibility  of  taking  more  corrector  iterations,  which  would  lead 
to  a  different  form  of  $tc  than  that  given  by  (1.2).  This  may  improve  stability  for 
some  problems  [16]. 

1.2.  Time  Step  Control.  We  assume  that  8n  is  computed  by  some  formula  like 
the  “switched  evolution  relaxation”  (SER)  method,  so  named  in  [20],  and  used  in, 
e.g.  [18],  [23],  and  [32].  In  its  simplest,  unprotected  form,  SER  increases  the  timestep 
in  inverse  proportion  to  the  residual  reduction. 

(1.5)  Sn  =  5ol|F(^o)||/||F(i„)||. 

(1.5)  implies  that,  for  n  >  1, 

_  \\FiXr.-l)\\ 

||F((r„)|i  • 

In  some  work  [16],  is  kept  below  a  large,  finite  bound  6max-  Sometimes  is  set 
to  c»  (called  “switchover  to  steady-state  form”  in  [13])  when  the  computed  value  of  8n 
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exceeds  8max-  In  view  of  these  practices,  we  will  allow  for  somewhat  more  generality 
in  the  formulation  of  the  sequence  {^„}.  We  will  assume  that  is  given  and  that 


(1.6) 


8n  =  4> 


(c 


for  n  >  1.  Here  ^  is  an  increasing  function  with  values  in  (0,  oo].  The  choice  in  [23] 
and  [32]  (equation  (1.5))  is 


m = e- 


Other  choices  could  either  limit  the  growth  of  S  or  allow  6  to  become  infinite  after 
finitely  many  steps.  One  variation  that  allows  for  all  of  these  possibilities  is 

So  if  =  oo  then  <f>{0  =  If  6  =  ^max/oi  <  oo  then 

?^(^)  —  niin(Q;^,  ^n^a^:) 


and  the  timesteps  are  held  bounded.  If  <  oo  and  S^ax  =  oo  then  switchover  to 
steady-state  form  is  permitted  after  a  finite  number  of  timesteps. 

In  [16]  the  timesteps  are  based  not  on  the  norms  of  the  nonlinear  residuals 
||F(a:„)||  but  on  the  norms  of  the  steps  — a;„_i||.  This  policy  has  risks  in  that  small 

steps  need  not  imply  small  residuals  or  accurate  results.  However  if  the  Jacobians 
are  uniformly  well  conditioned,  then  small  steps  do  imply  that  the  iteration  is  near  a 
root.  Here  formulae  of  the  type 

(1.8)  =  (?i>(^„_il|a:„  -  a;„_i||“^) 

are  used,  where  <j)  has  the  properties  mentioned  above.  The  function  given  by  (1.7), 
for  example,  is  a  reasonable  candidate. 

Our  formal  assumptions  on  the  function  <j)  are 
Assumption  1.1. 

1.  (f>  is  an  increasing  function  valued  in  (0,  ooj. 

2.  (j)  has  a  maximum  value,  S^ax  €  (0,  oo]  which  is  attained  at  €  (0,oo]. 

3.  >  a(  for  all  (  €  (O,^]- 

4.  Either  =  hmaxlo  or  <  00  and  Smax  =  00  • 

1.3.  Iteration  Phases.  We  divide  the  ^'tc  iteration  into  three  conceptually 
different  and  separately  addressed  phases. 

1.  The  initial  phase.  Here  6  is  small  and  x  is  far  from  a  steady  state  solution. 
This  phase  is  analyzed  in  §  2.3.  Success  in  this  phase  is  governed  by  stability 
and  accuracy  of  the  temporal  integration  scheme  and  proper  choice  of  initial 
data. 
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2.  The  midrange.  This  begins  with  an  accurate  solution  x  and  a  timestep  8 
that  may  be  small  and  produces  an  accurate  x  and  a  large  8.  We  analyze 
this  in  §  2.2.  To  allow  8  to  grow  without  loss  of  accuracy  in  x  we  make  a 
linear  stability  assumption  (part  4  of  Assumption  2.1). 

3.  The  terminal  phase.  Here  8  is  large  and  x  is  near  a  steady  state  solution. 
This  is  a  small  part  of  the  iteration,  usually  requiring  only  a  few  iterations. 
Aside  from  the  attention  that  must  be  paid  to  the  updating  rules  for  8^  the 
iteration  is  a  variation  of  the  chord  method  [24,  17]. 

We  analyze  the  terminal  phase  first,  as  it  is  the  easiest,  in  §  2.1.  Unlike  the  other 
two  phases,  the  analysis  of  the  terminal  phase  does  not  depend  on  the  dynamics 
of  x'  =  —V~^F{x).  The  initial  and  midrange  phases  are  considered  in  §  2.3,  with 
the  midrange  phase  considered  first  to  motivate  the  goals  of  the  initial  phase.  This 
decomposition  is  similar  to  that  proposed  for  GMRES  and  related  iterations  in  [21] 
and  is  supported  by  the  residual  norm  plots  reported  in  [23]  and  [10]. 

2.  Exact  Newton  Iteration.  In  this  section  we  analyze  the  three  phases  of 
the  solver  in  reverse  order.  This  ordering  makes  it  clear  how  the  output  of  an  earlier 
phase  must  conform  to  the  demands  of  the  later  phase. 

2.1.  Local  Convergence:  Terminal  Phase.  The  terminal  phase  of  the  iter¬ 
ation  can  be  analyzed  without  use  of  the  differential  equation  at  all. 

Lemma  2.1.  Let  {^„}  be  given  by  either  (1.6)  or  (1.8)  and  let  Assumption  1.1 
hold  with  a  <  1.  Let  F{x*)  =  0,  F\x*)  be  nonsingular,  and  F'  be  Lipschitz  continuous 
with  Lipschitz  constant  7  in  a  ball  of  radius  e  about  x*. 

Then  there  are  ei  >  0  and  Aq,  such  that  if  8max,8o  >  Aq  and  ||a:o  —  a;*||  <  Ci, 
then  the  sequence  defined  by  (1.2)  and  (1.6)  satisfies 

8n  ^  8fn,axt 

and  Xn  X*  q-superlinearly  if  8max  =  00  and  q-linearly  if  8max  <  00. 

Proof  Let  e  —  x  —  x*  denote  the  error.  As  is  standard,  [12],  [17],  we  analyze 
convergence  in  terms  of  the  transition  from  a  current  iterate  Xc  to  a  new  iterate  x+. 
We  must  also  keep  track  of  the  change  in  8  and  will  let  8c  and  6+  be  the  current  and 
new  pseudo-timesteps. 

The  standard  analysis  of  the  chord  method  ([17],  p.  76)  implies  that  there  are 
ei  <  e  and  A'c  such  that  if  ||ec||  <  ei, 

(2.1)  ||e+||</ic(|iec|P  +  ^;'||ee||). 

So  if  8c  >  Ao  and 


Cl  -t-  Ao  ^  <  ll{2Kc)i 
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then  ||e+||  <  ||ec||/2  and,  in  particular,  F{x+)  is  defined.  Reduce  ei  and  increase  Aq 
if  needed  so  that 


(2.2) 


Cl  +  Ao  ^  <  a/{8K{F'{x*)))  and  ej  <  3/o;, 


where  «  denotes  condition  number  and  a  is  from  part  3  of  Assumption  1.1.  Equations 
(2.1)  and  (2.2)  imply  that 


(2.3) 


||e+||<|M(ex  +  Ao-^)< 


(2.4) 


k{F'{x*)Y 

If  {6„}  is  computed  with  (1.6)  we  use  the  following  inequality  from  [17]  (p.  72) 

l|e+|| 


4Ac(F'(x*))||ec||  ||i^(a;c 


,  ||F(a:.,.)||  ^  4^(FV))||e^|| 


and  (2.3)  to  obtain 


> 


||F(:r+)||  -  4Ac(F'((r*))||e+| 
We  then  have  by  Assumption  1.1  that 


>  2a-\ 


8^  =  o6(6,||F((r,)||/||F(a:+)||)  >  mdc^)  >  <! 

where  is  from  part  2  of  Assumption  1.1. 

If  {^„)  is  computed  with  (1.8),  we  note  that 


^maxj  8i^Ct  ^ 
(  28c,  8ca  < 


|x+  -  Xell  <  lie+ll  +  lleell  <  3||ec||/2  <  3ej2  <  a/2 


and  hence 


=  (t>{8c/\\x+  -  Xc||)  >  <j){28cloi)  >  < 


^max^ 

28c, 


^cOC  >  6 

8cOC  < 


as  before. 

In  either  case,  >  dc  >  Aq  and  ||e+||  <  ||ec||/2.  Therefore  we  may  continue  the 
iteration  and  conclude  that  >  8max  and  ||e„||  0  at  least  q-linearly  with  q-factor 

of  1/2. 

If  8max  —  0O5  we  complete  the  proof  by  observing  that  since  — >  00  and  —>■  x* 

q-superlinear  convergence  follows  from  (2.1).  □ 

The  following  simple  corollary  of  (2.1)  applies  to  the  choice  =  a^. 

Corollary  2.2.  Let  the  assumptions  of  Lemma  2.1  hold.  Assume  that  = 
oc(-  ||eo||  <  Cl,  ^0  >  Ao,  and  8max  =  oo-  Then  the  convergence  of  {xn}  to  x*  is 
q-quadratic. 
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2.2.  Increasing  the  Time  Step:  Midrange  Phase.  Lemma  2.1  states  that 
the  second  phase  of  ^'tc  should  produce  both  a  large  S  and  an  accurate  solution.  We 
show  how  this  happens  if  the  initial  phase  can  provide  only  an  accurate  solution. 
This  clarifies  the  role  of  the  second  phase  in  increasing  the  timestep.  We  do  this  by 
showing  that  if  the  steady  state  problem  has  a  stable  solution,  then  'J’tc  behaves  well. 
We  now  make  assumptions  that  not  only  involve  the  nonlinear  function  but  also  the 
initial  data  and  the  dynamics  of  the  I  VP  (1.1). 

Assumption  2.1. 

1.  There  is  a  root-x*  of  F  that  satisfies  the  standard  assumptions  and  p  >  0 
such  that  if  \\z  —  xqH  <  V  then  the  solution  of  the  initial  value  problem 

(2.5)  x'  =  -V~^F{x),  x(0)  =  z 

converges  to  x*  as  t  oo. 

2.  F  is  everywhere  defined  and  Lipschitz  continuously  Frechet  differentiable 

3.  ||F'(x)||  <  M  for  all  x,  for  some  M  >  0 

4-  There  are  €2,^  >  0  such  that  if  \\x  —  x*||  <  €2  then  ||(/  +  6V“^F'(x))~^||  < 
(1  +  for  all  6  >0. 

The  analysis  of  the  midrange  uses  part  4  of  Assumption  2.1  in  an  important  way 
to  guarantee  stability.  The  method  for  updating  S  is  not  important  for  this  result. 

Theorem  2.3.  Let  {6„}  be  given  by  either  (1.6)  or  (1.8)  and  let  Assumption  1.1 
hold  with  a  <  1.  Let  Assumption  2.1  hold.  Let  6max  be  large  enough  for  the  conclu¬ 
sions  of  Lemma  2.1  to  hold.  Then  there  is  an  63  >  0  such  that  ^/||xo  — x*||  <  63,  and 
So  >  0,  either 


inf  Sn  =  0 

n 

or  Xji  ^  X  and  Sji  ^  Sj^mx. 

Proof.  Let  €3  <  min(ei,€2)  where  Ci  is  from  Lemma  2.1  and  62  is  from  part  4  of 
Assumption  2.1.  Note  that 

ei  =  eo-{I  +  6oV-^F'{xo))-^SoV-'^F{xo) 

=  eo  -  (/  +  SoV-^F'{xo))-^SoV-^F'{xo)eo 

+{I  +  6oV-^F'{xo))-^SoV-\F’{xo)eo  -  F{xo)) 

=  (/  +  6oV-^F'{xo))-^eo  +  (/  +  6oV-^F'{xo))-HoV-\F'{xo)eo  -  F{xo)) 
Now  there  is  a  c  >  0  such  that 

||V-^(F'(x)e-F(x))||<c||e||2 
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for  all  X  such  that  ||e||  <  ci-  Hence,  reducing  £3  further  if  needed  so  that  C3  < 
we  have 


(2.6) 


kill  <  Ikoll 


/l  +  ce3(5'\ 


<  Ikoll 


(i±^\ 

V  y 


<  Ikoll. 


If  inf„  ^„  =  ^*  >  0  then 


< 


'I +  ^8*  12' 
1  +  ^8* 


1 


for  all  n  >  1  and  hence  Xn  converges  to  x*  q-linearly  with  q-factor  {\+ ^8* l2){l-\- ^8*), 

This  convergence  implies  that  >  0  and  F{xn)  — >  0  so  (5n  — >■  8max  if  either  (1.6) 
or  (1.8)  is  used.  0 

This  result  says  that  once  the  iteration  has  found  a  sufficiently  good  solution, 
either  convergence  will  happen  or  the  the  iteration  will  stagnate  with  inf  =  0.  This 
latter  failure  mode  is,  of  course,  easy  to  detect.  Moreover,  the  radius  €3  of  the  ball 
about  the  root  in  Theorem  2.3  does  not  depend  on  inf 

2.3.  Integration  to  Steady  State:  Initial  Phase.  Theorem  2.3  requires  an 
accurate  estimate  of  x*,  but  asks  nothing  of  the  timestep.  In  this  section  we  show  that 
if  80  is  sufficiently  small,  and  (1.6)  is  used  to  update  the  timestep,  then  the  dynamics 
of  (1.1)  will  be  tracked  sufficiently  well  for  such  an  approximate  solution  to  be  found. 
It  is  not  clear  how  (1.8)  allows  for  this. 

Theorem  2.4.  Let  {6„}  6e  given  by  (1.6)  and  let  Assumption  1.1  hold  with 
o  <  1.  Let  Assumption  2.1  hold.  Let  e  >  0.  There  is  a  8  such  that  if  80  <8  then 
there  is  an  n  such  that  ||e„||  <  e. 

Proof.  Let  S  be  the  trajectory  of  the  solution  to  (2.5).  By  Assumption  2.1  x* 
satisfies  the  standard  assumptions  and  therefore  there  are  64  and  e/  such  that  if 


k  —  y\\  ^  some  y  £  S  and  ||jP(a;)||  <  ey 


then  lla;  —  x*||  <  €3,  which  will  suffice  for  the  conclusions  of  Theorem  2.3  to  hold.  Let 

M  =  sup  ||jF(a:)||. 

||x-y|l<£4  for  some  y  £  S 

We  will  show  that  if  sufficiently  small,  then  ||xfc  —  ar(tfc)||  <  64  until  ||i^(x/5)ll  < 
ej.  By  (1.6) 

(2.7)  Cl8o  =  <5o||F(a:o)||/M<4<^o||F(xo)|l/ey  =  Cu8o 

as  long  as  ||F(a;)||  >  e/  and  Xk  is  within  £4  of  the  trajectory  S.  Let  T  >  0  be  such 
that  if  is  the  solution  to  (2.5),  t>T.,  and  ||a:  —  z{t)\\  <  €4  then  ||jP(x)||  <  cy. 
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Let  ^  be  the  solution  of  (2.5).  Consider  the  approximate  integration  of  (2.5)  by 
(1.2).  Set 

1=0 

If  ||x„-2(^„)||  <  64  and  ||F(x„)||  >  tj  then  (2.7)  holds.  This  cannot  happen  if  >  T 
and  therefore  the  proof  will  be  complete  if  we  can  show  that 

ll^n  ■^(^n)ll  ^4 

until  the  first  n  such  that  ||F(x)||  <  ey,  which  will  happen  for  some 

(2.8)  n  <  ClT/6o: 

Note  that  (1.2)  may  be  written  as 

(2.9)  X„+a  =  X„  -  ^„y-'F(x„)  +  [/-(/  +  6r,V-^F'{Xn))-'^]Sr.V-'^F{Xn). 

There  is  an  mi  such  that  the  last  term  in  (2.9)  satisfies,  for  6n  sufficiently  small  and 
ll^n  •^(^«)ll  ^4? 

(2.10)  ||[/  -  (/  +  SnV-^F'{Xr,))-'^]SnV-^F{Xn)\\  <  niySl 

Let  En  =  \\xn  —  ^(^n)ll-  Then  we  have,  by  our  assumptions  on  F,  that  there  is 
an  Ml  >  0  such  that 

(2.11)  \\Xn  +  6nV~'^F{Xn)  “  a:(tn)  -  x'(t„)6„||  <  Mi6nEn. 

Finally,  there  is  an  m2  such  that  for  for  6n  sufficiently  small  and  ||xn  —  2(tn)||  <  ^4 

||a:(t„  +  6n)  -  x{tn)  -  x'(4)^„||  <  m2Sl. 

Setting  M2  =  mi  +  m2  we  have  for  all  n  >.l  (as  long  as  is  sufficiently  small 
and  \\xn  •^(^n)ll  ^  ^4) 

Fn  <  (1  +  Mi6n-l)En-l  +  M2^^_l. 

As  long  as  (2.7)  holds,  this  implies  that 

En  ^  (1  +  MiCu^o)En-l  +  M2C^8q. 

Consequently,  as  is  standard  [14,  15], 

6oM2C'[/[exp(nMiC[/6o)  -  1] 
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and  using  (2.8), 

(2.12) 

^  ^  ,  M2Cu[exp{CLMrCuT)  -  1] 

E^<  So 

So  if 

Ml' 

(2.1.3) 

^  ^^M2Cu[exp{CLMiCuT)  -  1] 

then  |(xn  -  z{tn)\\  <  ^4  for  all  n  until  F{xn)  <  £/  or  >  T.  This  completes  the  proof. 

□ 

The  problem  with  application  of  this  proof  to  the  update  given  by  (1.8)  is  that 
bounds  on  6  like  (2.7)  do  not  follow  from  the  update  formula. 

3.  Inexact  Newton  Iteration.  In  this  section  we  look  at  ^tc  as  implemented 
in  practice.  There  are  two  significant  differences  between  the  simple  version  in  §  2 
and  realistic  implementations: 

1.  The  Frechet  derivative  6~W  -\-F'{xn)  is  not  recomputed  with  every  timestep. 

2.  The  equation  for  the  Newton  step  is  solved  only  inexactly. 

Before  showing  how  the  results  in  §  2  are  affected  by  these  differences,  we  provide 
some  more  detail  and  motivation. 

Item  1  is  subtle.  If  one  is  solving  the  equation  for  the  Newton  step  with  a  direct 
method,  then  evaluation  and  factorization  of  the  .Jacobian  matrix  is  not  done  at  every 
timestep.  This  is  a  common  feature  of  many  ODE  and  DAE  codes,  [29,  25,  26,  3]. 
Jacobian  updating  is  an  issue  in  continuation  methods  [30,  27],  and  implementations 
of  the  chord  and  Shamanskii  [28]  methods  for  general  nonlinear  equations  [2,  17, 
24].  When  the  Jacobian  is  slowly  varying  as  a  function  of  time  or  the  continuation 
parameter,  sporadic  updating  of  the  Jacobian  leads  to  significant  performance  gains. 
One  must  decide  when  to  evaluate  and  factor  the  Jacobian  using  iteration  statistics 
and  (in  the  ODE  and  DAE  case)  estimates  of  truncation  error.  Temporal  truncation 
error  is  not  of  interest  to  us,  of  course,  if  we  seek  only  the  steady-state  solution. 

In  [16]  a  Jacobian  corresponding  to  a  lower-order  discretization  than  that  for  the 
residual  was  used  in  the  early  phases  of  the  iteration  and  in  [18],  in  the  context  of  a 
matrix-free  Newton  method,  the  same  was  used  as  a  preconditioner. 

'  The  risks  in  the  use  of  inaccurate  Jacobian  information  are  that  termination  deci¬ 
sions  for  the  Newton  iteration  and  the  decision  to  reevaluate  and  refactor  the  Jacobian 
are  related  and  one  can  be  misled  by  rapidly  varying  and  ill-conditioned  Jacobians 
into  premature  termination  of  the  nonlinear  iteration  [29,  31].  In  the  case  of  itera¬ 
tive  methods,  item  1  should  be  interpreted  to  mean  that  preconditioning  information 
(such  as  an  incomplete  factorization)  is  not  computed  at  every  timestep. 

Item  2  means  that  the  equation  for  the  Newton  step  is  solved  inexactly  in  the 
sense  of  [11],  so  that  instead  of 

(3.1) 


=  Xc  -t-  s 
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where  s  is  given  by  (1-3),  step  s  satisfies 


(3.2)  ||(r^F  +  F'ix))s  +  F(a;)||  <  »;||F(a:)||, 

for  some  small  ?/,  which  may  change  as  the  iteration  progresses.  Item  1  can  also  be 
viewed  as  an  inexact  Newton  method  with  rj  reflecting  the  difference  between  the 
approximate  and  exact  Jacobians. 

The  theory  in  §  2  is  not  changed  much  if  inexact  computation  of  the  step  is 
allowed.  The  proof  of  Lemma  2.1  is  affected  in  (2.1),  which  must  be  changed  to 

(3.3)  l|e+||</Cc(||e.ir  +  [i;‘+^J||e.||). 

This  changes  the  statement  of  the  lemma  to 

Lemma  3.1.  Let  {^„}  be  given  by  either  (1.6)  or  (1.8)  and  let  Assumption  1.1 
hold  with  a  <  1.  Let  F(x*)  =  0,  F'(x*)  be  nonsingular,  and  F'  be  Lipschitz  continuous 
with  Lipschitz  constant  7  in  a  ball  of  radius  e  about  x*. 

Then  there  are  ei>  0,  f}  and  Aq  such  that  if  Smaxi  bo  >  Aq,  rin'>fj  for  all  n,  and 
||a:o  —  ^*11  <  Iben  the  sequence  defined  by  (3.1),  (3.2),  and  (1.6)  satisfies 

bn  *■  bjnax 

and  Xn  —*■  X*  q-superlinearly  if  bmax  =  00  and  r/n  — ^  0  and  q-linearly  if  6max  <  00. 
Corollary  2.2  becomes 

Corollary  3.2.  Let  the  assumptions  of  Lemma  3.1  hold.  Assume  that  ^(if)  = 
af.  l|eo||  <  Cl,  So  >  Aq,  and  Smax  =  00.  Then  the  convergence  of  to  x*  is 
q-superlinear  ifrjn  0,  and  locally  q-quadratic  if  qn  =  0(||F(x„)||). 

The  analysis  of  the  midrange  phase  changes  in  (2.6),  where  we  obtain 

(3.4)  ||ei||  <  ||eo||  (^Kr,Vo  + 

for  some  /i'„  >  0.  This  means  that  fj  must  be  small  enough  to  maintain  the  q-linear 
convergence  of  {x„}  during  this  phase.  The  inexact  form  of  Theorem  2.3  is 

Theorem  3.3.  Let  {x„}  be  given  by  (3.1)  and  (3.2)  and  let  be  given  by 
either  (1.6)  or  (1.8).  Let  Assumption  1.1  hold  with  a  <  1.  Let  Assumption  2.1  hold. 
Let  bmax  be  large  enough  for  the  conclusions  of  Lemma  2.1  to  hold.  Then  there  are 
€3  >  0  and  fj  such  that  if  Pn  <  fj,  ||a:o  -  ic*||  <  and  So  >  0,  either 

inf  =  0 

n 

or  Xn  ^  X  and  Sn  ^  Smax* 

Inexact  Newton  methods,  in  particular  Newton- Krylov  solvers,  have  been  applied 
to  ODE/DAE  solvers  in  [1],  [5],  [4],  [6],  [7],  and  [9].  The  context  here  is  different  in 
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that  the  nonlinear  residual  F{x)  does  not  reflect  the  error  in  the  transient  solution 
but  in  the  steady  state  solution. 

The  analysis  of  the  initial  phase  changes  through  (2.10).  We  must  now  estimate 

6nV~'^F{Xn)  +  5n- 

Set  r  =  +  F'{x.n))sn  +  F(xn).  Note  that 

=6ry-\F{Xr,)  +  8-^VSn) 

=  8ry-\r  -  F'{Xn)Sn). 

Now, 

=  {I+8nV-^F\Xr:))-Hry-\r  -  F{Xn)) 

and  hence,  assuming  that  the  operators  {I F  8nV~'^F'{xn))~'^  are  uniformly  bounded, 
there  is  m3  such  that 


and  hence 

(3.5)  \\8nV-^Fixr.)  +  3„|1  <  ^„l|V-'||(r;„||F(x„)||  +  \\F'{xr,)\\Tm8n). 

We  express  (3.5)  as 

(3.6)  \\8nV-'^F{xn)  +  5„||  <  mi{rin8n  +  8l). 

Hence,  if  =  O(8o)  or  rjn  =  0{En)  we  still  obtain  (2.12)  and  the  inexact  form 
of  Theorem  2.4: 

Theorem  3.4.  Let  {x„}  be  given  by  (3.1)  and  (3.2)  and  let  {^„}  be  given  by 

(1.6) .  Let  Assumption  1.1  hold  with  a  <  1.  Let  Assumption  2.1  hold.  Assume  that 
the  operators  (/  +  8nV~^F'{xn))~^  are  uniformly  bounded  in  n.  Let  e>  0.  There  are 
8  and  fj  such  that  if  8o  <  8  and  ?7n  <  ^  then  there  is  an  n  such  that  ||e„||  <  e. 

The  restrictions  on  rj  in  Theorem  3.4  seem  to  be  stronger  than  those  on  the  results 
on  the  midrange  and  terminal  phases.  This  is  consistent  with  the  tight  defaults  on  the 
forcing  terms  for  Newton-Krylov  methods  when  applied  in  the  ODE/DAE  context 
[1,  5,  6,  7,  9]. 

4.  Numerical  Experiments.  In  this  section  we  examine  a  commonly  used  ^tc 
technique,  switched  evolution/relaxation  (SER)  [20],  applied  to  a  Newton-like  method 
for  inviscid  compressible  flow  over  a  four-element  airfoil  in  two  dimensions.  Three 
phases  corresponding  roughly  to  the  theoretically-motivated  iteration  phases  of  §  2 
may  be  identified.  We  also  compare  SER  with  a  different  ^tc  technique  based  on 
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Fig.  4.1.  Unstructured  grid  around  four- element  airfoil  in  landing  configuration  —  near-field 
view, 

bounding  temporal  truncation  error  (TTE)  [19].  TTE  is  slightly  more  aggressive 
than  SER  in  building  up  the  time  step  in  this  particular  problem,  but  the  behavior 
of  the  system  is  qualitatively  the  same. 

The  physical  problem,  its  discretization,  and  its  algorithmic  treatment  in  both  a 
nonlinear  defect  correction  iteration  and  in  a  Newton-Krylov-Schwarz  iteration  —  as 
well  as  its  suitability  for  parallel  computation  —  have  been  documented  in  earlier  pa¬ 
pers,  most  recently  [10]  and  the  references  therein.  Our  description  is  correspondingly 
compact. 

The  unknowns  of  the  problem  are  nodal  values  of  the  fluid  density,  velocities, 
and  specific  total  energy,  x  =  (. . . ,  pi,  Ui,  v,,  e,-, . .  .)^  at  N  vertices  in  an  unstructured 
grid  of  triangular  cells  (see  Fig.  4.1).  The  system  F{x)  =  0  is  a  discretization  of  the 
steady  Euler  equations: 


(4.1) 

V  •  (pv)  =  0 

(4,2) 

V  •  (pvv  +  pl)  =  0 

(4.3) 

V  •  ((pe -H  p)v)  =  0 

12 

where  the  pressure  p  is  supplied  from  the  ideal  gas  law,  p  =  ^0(7-  l)(e-  |vp/2),  and  7 
is  the  ratio  of  specific  heats.  The  discretization  is  based  on  a  control  volume  approach, 
in  which  the  control  volumes  are  the  duals  of  the  triangular  cells  —  nonoverlapping 
polygons  surrounding  each  vertex  whose  perimeter  segments  join  adjacent  cell  centers 
to  midpoints  of  incident  cell  edges.  Integrals  of  (4.1)-(4.3)  over  the  control  volumes 
are  transformed  by  the  divergence  theorem  to  contour  integrals  of  fluxes,  which  are 
estimated  numerically  through  an  upwind  scheme  of  Roe  type.  The  effective  scaling 
matrix  V  for  the  ^tc  term  is  a  diagonal  matrix  that  depends  upon  the  mesh. 

The  boundary  conditions  correspond  to  landing  configuration  conditions:  sub¬ 
sonic  (Mach  number  of  0.2)  with  a  high  angle  of  attack  of  (5®).  The  full  adaptively 
clustered  unstructured  grid  contains  6,019  vertices,  with  four  degrees  of  freedom  per 
vertex  (giving  24,076  as  the  algebraic  dimension  of  the  discrete  nonlinear  problem). 
Figure  4.1  shows  only  a  near-held  zoom  on  the  full  grid,  whose  far-held  boundaries 
are  approximately  twenty  chords  away.  The  initial  pseudo-timestep,  60  «  3  x  10“^ 
corresponds  to  a  CFL  number  of  20.  The  pseudo-timestep  is  allowed  to  grow  up  to 
six  orders  of  magnitude  over  the  course  of  the  iterations.  It  is  ultimately  bounded  at 
^max  =  10®  ■  ^0  guaranteeing  a  modest  diagonal  contribution  that  aids  the  invertibility 
of  {6-^V  +  F'{Xr,))-K 

The  initial  iterate  is  a  uniform  how,  based  on  the  far  held  boundary  conditions 
—  constant  density  and  energy,  and  constant  velocity  at  a  given  angle  of  attack. 

The  solution  algorithm  is  a  hybrid  between  a  nonlinear  defect  correction  and  a  full 
Newton  method,  a  distinction  which  requires  further  discussion  of  the  processes  that 
supply  F{x)  and  F'{x)  within  the  code.  The  form  of  the  vector-valued  function  F{x) 
determines  the  quality  of  the  solution  and  is  always  discretized  to  required  accuracy 
(second-order  in  this  paper).  The  form  of  the  approximate  Jacobian  matrix  F'{x), 
together  with  the  scaling  matrix  V  and  time  step  6^  determines  the  rate  at  which 
the  solution  is  achieved  but  does  not  affect  the  quality  of  a  converged  result,  and 
is  therefore  a  candidate  for  replacement  with  a  matrix  that  is  more  convenient.  In 
practice,  we  perform  the  matrix  inversion  in  (1.2)  by  Krylov  iteration,  which  requires 
only  the  action  of  F'{x)  on  a  series  of  Krylov  vectors,  and  not  the  actual  elements 
of  F'{x).  The  Krylov  method  was  restarted  GMRES(20)  preconditioned  with  1-cell 
overlap  additive  Schwarz  (8  subdomains). 

Following  [5, 8],  we  use  matrix-free  Frechet  approximations  of  the  required  action: 

(4.4)  F'{x)v  ~  —  [F{x  -t-  hv)  —  F(a:)] . 

However,  when  preconditioning  the  solution  of  (1.2),  we  use  a  more  economical  matrix 
than  the  Jacobian  based  on  the  true  F{x),  obtained  from  a  first-order  discretization 
of  the  governing  Euler  system.  This  not  only  decreases  the  number  of  elements  in 
the  preconditioner,  relative  to  a  true  Jacobian,  but  also  the  computation  and  (in 
the  parallel  context)  communication  in  applying  the  preconditioner.  It  also  results 
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Residual,  Update,  and  Timestep  for  SER 


Fig.  4.2.  SER  convergence  history 

in  a  more  numerically  diffusive  and  stable  matrix,  which  is  desirable  for  inversion. 
The  price  for  these  advantages  is  that  the  preconditioning  is  inconsistent  with  the 
true  Jacobian,  so  more  linear  subiterations  may  be  required  to  meet  a  given  linear 
convergence  tolerance.  This  has  an  indirect  feedback  on  the  nonlinear  convergence 
rate,  since  we  limit  the  work  performed  in  any  linear  subiteration  in  an  inexact  Newton 
sense. 

In  previous  work  on  the  Euler  and  Navier-Stokes  equations  [10,  22],  we  have  noted 
that  a  ^tc  method  based  on  a  consistent  high-order  Jacobian  stumbles  around  with  a 
nonmonotonic  steady-state  residual  norm  at  the  outset  of  the  nonlinear  iterations  for 
a  typical  convenient  initial  iterate  far  from  the  steady-state  solution.  On  the  other 
hand,  a  simple  defect  correction  approach,  in  which  F'{x)  is  based  on  a  regularizing 
first-order  discretization  everywhere  it  appears  in  the  solution  of  (1.2),  not  just  in  the 
preconditioning,  allows  the  residual  to  drop  smoothly  from  the  outset.  In  this  work, 
we  employ  a  hybrid  strategy,  in  which  defect  correction  is  used  until  the  residual  norm 
has  fallen  by  three  orders  of  magnitude,  and  inexact  Newton  thereafter.  As  noted  in 
§  .3,  inexact  iteration  based  on  the  true  Jacobian  and  iteration  with  an  inconsistent 
Jacobian  can  both  be  gathered  under  the  rj  of  (3.2),  so  the  theory  extends  in  principal 
to  both. 

With  this  background  we  turn  our  attention  to  Fig.  4.2,  in  which  are  plotted 
on  a  logarithmic  scale  against  the  $tc  iteration  number:  the  steady-state  residual 
norm  ||F(a:„)||2  at  the  beginning  of  each  iteration,  the  norm  of  the  update  vector 
l|x„+i  —  Xn\\2,  and  the  pseudo-timestep  6n- 

The  residual  norm  falls  nearly  monotonically,  as  does  the  norm  of  the  solution 
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update.  Asymptotic  convergence  cannot  be  expected  to  be  quadratic  or  superlinear, 
since  we  do  not  enforce  7/„  — ^  0  in  (3.5).  However,  linear  convergence  is  steep,  and  our 
experience  shows  that  overall  execution  time  is  increased  if  too  many  linear  iterations 
are  employed  in  order  to  enforce  0  asymptotically.  In  the  results  shown  in  this 
section,  the  inner  linear  convergence  tolerance  was  set  at  10“^  for  the  defect  correction 
part  of  the  trajectory,  and  at  10"^  for  the  Newton  part.  The  work  was  also  limited 
to  a  maximum  of  12  restart  cycles  of  20  Krylov  vectors  each. 

Examination  of  the  pseudo-timestep  history  shows  monotonic  growth  that  is  grad¬ 
ual  through  the  defect  correction  phase  (ending  at  n  =  14),  then  more  rapidly  grow¬ 
ing,  and  asymptotically  at  dmax  (beginning  at  n  =  20).  Steps  21,  24,  and  32  show 
momentary  retreats  from  8max  in  response  to  a  refinement  on  the  $tc  strategy  that 
automatically  cuts  back  the  pseudo-timestep  by  a  fixed  factor  if  a  nonlinear  resid¬ 
ual  reduction  of  less  than  |  is  achieved  at  the  exhaustion  of  the  maximum  number 
of  Krylov  restarts  in  the  previous  step  (during  the  terminal  Newton  phase).  Close 
examination  usually  reveals  a  stagnation  plateau  in  the  linear  solver,  and  it  is  more 
cost  effective  to  fall  back  to  the  physical  transient  to  proceed  than  to  grind  on  the 
ill-conditioned  linear  problem.  These  glitches  in  the  convergence  of  ||F(a;„)||2  are  not 
of  nonlinear  origin. 

Another  timestep  policy,  common  in  the  ODE  literature,  is  based  on  controlling 
temporal  truncation  error  estimates.  Though  we  do  not  need  to  maintain  temporal 
truncation  errors  at  low  levels  when  we  are  not  attempting  to  follow  physical  tran¬ 
sients,  we  may  maintain  them  at  high  levels  as  a  heuristic  form  of  stepsize  control. 
This  policy  seems  rare  in  external  aerodynamic  simulations,  but  is  popular  in  the 
combustion  community  and  is  implemented  in  [19].  The  first  neglected  term  in  the 
Euler  discretization  of  ^  is  \x"  •  so  a  reasonable  mixed  absolute- relative  bound 
on  the  error  in  the  component  of  x  at  the  step  is 


(4.5) 


81  d-^x, 

2(1  +  la:,])  dt^ 


<  T, 


where  (a:f)n  can  be  approximated  by 


2 

8n-l  +  8n-2 


{^i)n  {^i)n-l 

8n-l 


(Xj)n— 1  (^i)n— 2 

8n-2 


Taking  r  as  |  and  implementing  this  strategy  in  the  Euler  code  in  place  of  SER 
yields  the  results  in  Fig.  4.3.  Arrival  at  8max  occurs  at  the  same  step  as  for  SER,  and 
arrival  at  the  threshold  ||F’(a;„)||  <  10"^^  occurs  one  iteration  earlier.  However,  the 
convergence  difficulties  after  having  arrived  at  8max  are  slightly  greater. 

5.  Conclusions.  Though  the  numerical  experiments  of  the  previous  section  do 
not  confirm  the  theory  in  detail,  in  the  sense  that  we  do  not  verify  the  estimates  in  the 
hypotheses,  a  reassuring  similarity  exists  between  the  observations  of  the  numerics 
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Residual,  Update,  and  Timestep  for  TTE 


Fig.  4.3.  TTE  convergence  history 


and  the  conceptual  framework  of  the  theory,  which  was  originally  motivated  by  similar 
observations  in  the  literature.  There  is  a  fairly  long  induction  phase,  in  which  the 
initial  iterate  is  guided  towards  the  Newton  convergence  domain  by  remaining  close 
to  the  physical  transient,  with  relatively  small  timesteps.  There  is  a  terminal  phase 
which  can  be  made  as  rapid  as  the  capability  of  the  linear  solver  permits  (which 
varies  from  application  to  application),  in  which  an  iterate  in  the  Newton  convergence 
domain  is  polished.  Connecting  the  two  is  a  phase  of  moderate  length  during  which 
the  time  step  is  built  up  towards  the  Newton  limit  of  Smax,  starting  from  a  reasonably 
accurate  iterate.  The  division  between  these  phases  is  not  always  clear  cut,  though 
exogenous  experience  suggests  that  it  becomes  more  so  when  the  corrector  of  §  1  is 
iterated  towards  convergence  on  each  time  step.  We  plan  to  examine  this  region  of 
parameter  space  in  conjunction  with  an  extension  of  the  theory  to  mixed  steady/^tc 
systems  (analogs  of  differential-algebraic  systems  in  the  ODE  context)  in  the  future. 
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